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A TAYLOR-GALERKIN FINITE ELEMENT ALGORITHM FOR TRANSIENT 
NONLINEAR THERMAL-STRUCTURAL ANALYSIS 


Earl A. Thornton and Pramote Dechaumphai** 
Old Dominion University 
Norfolk, Virginia 


Abstract 


A Taylor-Galerkin finite element method for 
solving large, nonlinear thermal -structural 
problems is presented. The algorithm is 
formulated for coupled transient and uncoupled 
quasistatic thermal -structural problems. 
Vectorizing strategies ensure computational 
efficiency. Two applications demonstrate the 
validity of the approach for analyzing transient 
and quasistatic thermal -structural problems. 
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displacement components 
velocity components 
typical unknown 
coordinate directions 

thermal-mechanical coupling 
coefficients, eq. (lc) 

densi ty 
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Nomenclature 

propagation speed, eq. (16) 
element area 

fictitious damping constant, eq. 
(3a) 

elastic constants 

specific heat at constant volume 

diffusion coefficient, eq, (16) 

body forces per unit volume 

stresses or heat fluxes 

"load" vector, eq. (la) 

thermal conductivities 

components of unit normal surface 
vector 

mass matrix 

element interpolation functions 

interpolation functions on an 
element edge 

heat flux components 
volumetric heating rate 
load vectors, eq. (14) 
time 

time step 
temperature 

reference temperature for zero 
stress 
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Subscripts 

D constant element quantity 
S constant surface quantity 
n time step index 


Superscript 

n time step index 


Introduction 

The determination of the structural 
response induced by thermal effects is an 
important factor in many aerospace structural 
designs. Extreme aerodynamic heating on 
advanced aerospace vehicles may produce severe 
thermal stresses that can reduce operational 
performance or even damage structures. The 
performance of laser devices can be degraded by 
thermal distortions of mirror surfaces. The 
thermal environment in space may cause orbiting 
structures to distort beyond operational 
tolerances. To predict the structural response 
accurately, effective numerical techniques 
capable of both thermal and structural analyses 
are required. The finite element method has 
been found to be particularly suited for such 
analyses due to its capability to model complex 
geometry and to perform both thermal and 
structural analyses. 

Algorithms required for analyzing 
thermally-induced structural behavior depend on 
the rate at which structural temperatures vary 
with time. When temperatures change rapidly the 
analyses are strongly coupled, and thermal- 
structural interactions occur that mandate 
simultaneous solution of the thermal-structural 
equations. Temperature changes can occur 
rapidly due to propagation of thermoelastic 
waves, during vibration induced by periodic 
variations of temperature fields, due to thermal 
shocks and in similar circumstances. These 
types of problems may involve resolving wave- 
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like details of the time dependent response for 
complex structures. Moreover, if the mechanical 
deformation terms in the heat transfer energy 
equation are retained, the equations are 
inherently non-linear even in the material's 
elastic range. There is a need for effective 
finite • element solution algorithms that can 
solve large nonlinear transient problems 
efficiently using new vector or parallel 
computing technology. 

In the most common approach to determining 
structural responses Induced by thermal effects, 
the thermal-structural analyses are assumed 
uncoupled, and the structural analysis is 
assumed quasi static. The uncoupled assumption 
means that mechanical deformation terms in the 
heat transfer energy equation are neglected. 
The quasistatic assumption means that inertia 
terms in the structural equations of motion are 
neglected. The practical effect of these 

assumptions is that the heat transfer analysis 
can be performed first, and the resulting 
temperatures can be used as input to a 

subsequent stress analysis. This approach works 
well when temperatures change slowly as occurs, 
for example, in structures subjected to 
aerodynamic heating. Under these circumstances, 
the uncoupled, quasistatic idealization provides 

an effective approach for finite element 

thermal -structural analysis. 

The purpose of this paper is to present a 
Taylor-Galerkin finite element method for 
solving large, nonlinear transient thermal- 
structural problems. The method is an 
application of a Taylor-Galerkin algorithm 
recently developed to solve the conservation 
equations of inviscid, compressible flow 1-2 . In 
the flow problem, the algorithm is used to solve 
the highly nonlinear Euler equations that 

includes capturing shock discontinuities in the 
flow field. Finite element models of flow 

problems usually are quite large with the number 
of equations typically in the range from 3,000 
to 100,000 or more. Thus, the algorithm appears 
to have desirable attributes that will make it 
effective for large nonlinear, transient 
thermal -structural problems. 

The formulation of nonlinear thermal- 
structural problems will be presented first, 
then the Taylor-Galerkin algorithm will be 
described. Next, explicit evaluation of the 
finite element integrals is described. Then the 
programming strategy for a vector computer 
implementation of the algorithm is described. 
Finally, results from two thermal-structural 
applications are presented. 


Thermal-Structural Formulations 


Coupled Thermal -Structural Analyses 

The nonlinear coupled thermal-structural 
equations for a two-dimensional continuum 3 can 
be written in the form 


MU) jME) b{F } = 

5t bx by 


{H} 


' (la) 


where (U) is the vector of unknowns, (E) and 
(F) are vectors of stresses and heat fluxes, and 
{H > is a "load" vector. For coupled thermal- 
structural analyses a classical form of coupling 
is assumed in the equation for illustrative 
purposes, although more general thermal- 
structural coupling is permitted. 
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where u $nd v are the displacements, p is the 
density, u and v are the velocities, c v is the 
specific heat and T is the temperature; 
a x‘ a y’ T xy are tt1e stress components; 

q , q are the heat fluxes; f , f are body 

^x y x y 

force components per unit volume; p^ ^>2 

are coefficients that depend on the coefficients 
of thermal expansion, and Q is the internal heat 
generation rate per unit volume. The first two 
equations in eq. (1) define the velocity 
components, the third and fourth equations are 
the equations of motion, and the fifth equation 
represents conservation of energy. The term 
with the square brackets in the last line of 
(H) represents the classical nonlinear form of 
thermal-structural coupling 3 . 

The structural and thermal equations are 
written in the form of eq. (1) to resemble the 
conservation equations of fluid flow. In the 
fluid context, the components of (U) are called 
the conservation variables, and the components 
of (E) and (F) are fluxes of mass, momentum and 
energy across the faces of a control volume. In 
the thermal -structural context, the components 
of U may also be regarded as conservation 
variables. The stress components of (E) and 
(F> now represent tractions on surfaces of the 
control volume; however, q x and q y still 
represent heat fluxes across surfaces of the 
control volume. 

In formulating the constitutive equations, 
highly nonlinear relations between stresses. 
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strains and temperatures are permitted as well 
as nonlinear relations between heat fluxes, 
temperatures and temperature gradients. 
Anisotropic materials can be accommodated as 
well. For simplicity, simple constitutive 
relations for linear, elastic orthotropic 
materials will be illustrated. The stress- 
strain relations for an orthotropic material are 
expressed as 
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where C^j are elastic constants, and T 0 is the 
reference temperature for zero stress. The heat 
fluxes are expressed by Fourier's law. 
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pf x 

pf 


(3b) 


The stress and heat flux components have been 
defined in eq. (2) . 


Boundary and Initial Conditions 


Equations (1) and (3) are solved subject to 
appropriate initial and boundary conditions. 
The initial conditions consist of specifying the 
distributions of the conservation variables 
(U) at time zero. The structural boundary 
conditions consist of specifying the 
displacements or surface tractions at all points 
on the boundary. The thermal boundary 
conditions consist of specifying temperatures or 
heat fluxes at all points on the boundary. 
Convective and radiation boundary conditions are 
incorporated through heat fluxes. 
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q x = ' k x S 


k *1 

y ay 


(2b) 


where k x , and k are the thermal conductivities. 


The solution domain D is divided into an 
arbitrary number of elements of r nodes each. 
Figure 1 shows the thermal and structural 
models. The same finite element model is used 
for the thermal and structural analyses. The 
figure shows typical quadrilateral elements 
(r=4) used in this paper. For simplicity, the 
finite element formulation will be given for a 
single scalar equation. 


Quasistatic Thermal-Structural Analysis 
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When temperatures change slowly, the 
inertia terms in the structural equations of 
motion can be neglected. This means that the 

first two equations in Eq. (1) are not 

required. However, even for the quasistatic 
case, the equations will be solved using a time- 
marching scheme. To retain the thermal- 
structural equations in a form suitable for 
time-marching, the components of eq. (la) are 
written as 
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where c is a fictitious damping constant that is 
used to facilitate time marching to a steady- 
state quasistatic solution. The new load vector 
is 


where the variables u, E, F and H are analogous 
to the corresponding vector quantities in eq. 
(1). Let (u) n denote the element nodal values 
of the variable u(x,y,t) at time t n . The time 
step At spans two typical times t n and t n+1 in 
the transient response. The computation 
proceeds through two time levels, t n+1/ , 2 and 
t n+ i. At time level t n+ i/ 2 , vaiues for u that 
are constant within each element are computed 
explicitly. At time level t n+1 , the constant 
element values computed at the first time level 
are used to compute nodal values for u. In the 
time level tp^j computations, element contribu- 
tions are assembled to yield the global 
equations for nodal unknowns. The resulting 
equations are approximately diagonalized to 
yield an explicit algorithm. 

Time Level t n+1/2 


To advance the solution to time level 
t n+1 / 2 , a truncated Taylor series yields 
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u(x,y,t n+l/2 ) = u(x,y,t n ) (5) 

+ i At S (x * y *V* 

Then eq. (4) is introduced on the right hand 
side of eq. (5) so that 

u(x,y,t n+1/2 ) = u(x,y,t n ) 


it 


H (x » y *V + § (x * y,t n ) 


+ 2 At H(x,y.t n ) (6) 


With eq (8), the dependent variable for 

each element can be computed explicitly using 
nodal values for u, E, F and H from the previous 
time t n . A later section will discuss the 
evaluation of the three integrals that appear on 
the right hand side of eq. (8). 

In advancing the solution to the next time 
level, the values of the dependent variables on 
surfaces with specified flux boundary conditions 
may be required also. Let u s n+ ' 2 denote the 
surface value on a typical element edge IJ on 
the surface. Fig. 1. Following the approach 
used previously, u s n+ ^ 2 is assumed constant on 
edge IJ at time t n+J ^ 2 , but at time t n , u, E, F 
and H vary along the edge. Thus the following 
approximations are used on an element edge 


At time level t n+ i/ 2 * the dependent variable 
u(x,y,t n+1 ^ 2 ) is assumed to have a constant 

value Up 0 ^ 2 within an element. 

At time level i^, in the response u, E, F 
and H vary within an element and are 
interpolated from nodal values. Thus, the 
following spatial approximations are used within 


an element. 



u( x .y,t n+ i^ 2 ) 

. m n+1/2 
" u 0 

(7a) 

u(x,y,t n ) 

= [N(x.y)](u} n 

(7b) 

E(x,y,t n ) 

= [N(x,y)]{E) n 

(7c) 

F(x,y.t n ) 

• [N(x,y)](F} n 

(7d) 

H(x,y,t n ) 

= CN(x,y)]{H} n 

(7e) 


where [N(x,y>] denotes element interpolation 
functions and {u} n is a vector of the element 
nodal quantities. The equations for u 0 n+1//2 
for each element are derived by the method of 
weighted residual s A . The spatial approximations 
given in eq. (7) are introduced into eq. (6) to 
give a residual; the result is multiplied by a 
weighting function which in this case is 
unity. Finally, the weighted residual is 
integrated over the area A of the element. The 
result is, 

A U D n+1/2 = / A [N] dA {u}" 


r rSN 1 j * / r \ n At r r SN-. , . r j - a n 
T A c d7 ] dA {E} - — >a [ a7 ] dA {F} 


/ \ n+1/2 

u x, y ,t n+l/2 = u s 

(9a) 

u(x,y,t n ) = [N( s) ] (u } n 

(9b) 

E(x,y,t n ) = [N( s) ] (E ) n 

(9c ; 

F(x,y,t n ) = [N(s)](F} n 

(9d) 

H(x,y,t n ) = [N( s) ] (H } n 

( 9e ) 

where [N(s)] denotes the 

functions. Using the method 

residuals, the values for u s n+ ' 
integrating along the length L 

interpolation 
of weighted 
are derived by 
of an element 


edge. Hence 

n+1/2 , n At , 8N n 

L U S = CN3 ds (u } - — J L [— ] ds IE} 

- f L [|£] ds (F } n + f L [N] ds{H} n . (10) 

Thus, eqs. (8) and (10) can be used to advance 
explicitly the element and surface values of the 
dependent variables to t n+ j/ 2 . Beginning with 

nodal values of (u) n , (E} n , (F} n and (H) n , eq. 
(8) is used to compute constant values 

u D n+1/2 for each element. In a similar way, eq. 
(10) is used to compute constant surface values 
u n+ ^ for element boundary edges that require 
these values at t n+1 / 2 . These values are 

computed explicitly by looping through all 
elements and appropriate element edges. 

Time Level t n+ ^ 


+ ■j— [N] dA { H } . 


( 8 ) 


To advance the solution to t n+ j, forward 
and backward truncated Taylor series expansions 
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(15c) 


at t n+ */2 are used to write the approximation 


{R 2 ) n = - At / s (N }[N] ds (A {E}£ 


ujx.y.t^i) = u(x.y.t n ) (11) 

+ At at ^ x,y,t n+l/2^ ‘ 

Then, following the approach used previously, 
eq. (4) is introduced on the right side to yield 

u(x,y,t n+1 ) = u(x,y,t n ) (12) 
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The equation for the nodal values of (u) n+1 can 
next be derived by the method of weighted 
residuals using the interpolation functions 
as weighting functions. These terms containing 
the derivatives are Integrated by parts. The 
resulting equation contains values of the fluxes 
E and F evaluated at t n+1 / 2 * The fluxes at the 
mid-step are linearly interpolated from their 
values at t„ and t^. Thus 

E(x,y,t n+1 y 2 ) = (1-9) E(x,y,t n ) (13) 

+ 9 E(x » y -W 

F(x,y,t n+l/2 ) = (1 ' 9) Flx >y*V + 9 F ^ x »y» t n +i^ 

where the interpolation parameter 9 varies from 
zero to one. These operations yield the 
equations for the nodal values of a single 
element, 

[M]{u} n+1 = [M]{u} n + (1-0) {R 1 > n * (1-0) (R 2 > n 
+ 9{R 3 ) n+1 + e(R 4 ) n+1 + {R 5 ) n+1/2 (14) 


where 


[M] = J A (N)[N] dA (15a) 


{R } n = At /. A [N] dA (E) n 
1 A Oa 


(15b) 


In eqs. (15c) and (15e) the coefficients A and m 
are the components of a unit vector normal to 
the boundary, see Fig. 1. Following usual 
finite element procedures, the element matrices 
given in eq. (15) are assembled to form system 
equations. 


The matrix [M] defined by eq. (15a) is the 
element consistent mass matrix. The terms 
{R. }, i = 1,2, ...5 , defined by eqs. (15b - 15f) 
represent “load" vectors. The load vectors 
(R 3 > n and {R 2 } n are computed from value known 
at time step t^,; the load vector {Rg} n+li ' 2 is 


computed 

t n+l/2 • 
n+1 


from values 
The load 


computed at 
vectors 


time step 
(Rg } n+ ^ and 


(R 4 > 
at t, 


n+l‘ 


depend on values yet to be determined 
The consistent mass matrix and the 
latter two load vectors mean that the algorithm 
in the general form of eq. (14) with 

arbitrary 9 is an implicit scheme. For the 
computations presented in this paper, the 
algorithm was cast into an explicit form by 
using a lumped mass matrix and selecting 

9 = 0 to eliminate (R 3 ) n+ ^ and (R 4 ) n+ * . 


The explicit two time-level Taylor-Galerkin 
algorithm is conditionally stable. Guidelines 
for determining allowable time steps have been 
established by studying linear one-dimensional 
model equations selected to represent the 
mathematical character of the physical 
problems. An analysis of the coupled thermo- 
elasticity problem shows that the partial 
differential equations are hyperbolic. The 
stability of the algorithm was studied by 
considering the linear equation. 
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A 3y 


where a is a propagation 

speed, and D is a diffusion coefficient. A 
stability analysis of the algorithm applied to 
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this equation shows that the time step should 
sati sfy 


At < 


AX 
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1 + 2 
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aAx 

T3 


(17) 


where Ax is the mesh spacing. An analysis of 
the quasistatic thermal -structural problem shows 
that the partial differential equations are 
parabolic. The stability of the algorithm was 
studied by considering the linear parabolic 
equation. 


du 

at 


a 2 u 

17 


(18) 


A stability analysis of the algorithm applied to 
this equation shows that the time step should 
satisfy 


At < 


ax 2 

nr 


(19) 


For practical computations, eqs. (17) and (19) 
have been used as a guide to estimate the time 
step, but insight into the physical problem and 
some trial and error have been needed to perform 
the computations successfully. 


vectorization scheme are the number of elements 
in the finite element model with occasional 
operations using vector lengths equal to the 
number of nodes. 

The critical vectorization tasks are for 
operations that are repetitive and performed at 
every time step. For finite element algorithms, 
these operations are: (1) assembly of element 

contributions into the global system of 
equations, (2) solution of the global system of 
equations, and (3) application of boundary 
conditions. 

The assembly of element contributions into 
the global system of equations requires special 
routines for vectorization. Nodal unknowns are 
stored in one dimensional arrays from I' to the 
number of nodes in the model, and in general, 
node numbering may be arbitrary throughout the 
mesh. Assembly of element contributions is 
performed using the VPS 32 FORTRAN-suppl ied 
scatter routine which places an element 
contribution into the proper location in the 
system of equations based on the element 
connectivity. Every element that contains a 
particular node in its connectivity provides its 
own contribution to the system equations, 
therefore, the assembly is an additive 
operation. "Scattering" alone would merely 
overwrite a previous element contribution. The 
special vector routine, then does an "additive 
sea tter." 


Explicit Evaluation of Element Integrals 

Element integrals for the Taylor-Galerkin 
algorithm shown in eqs. 8, 10, and 15 were 
evaluated in closed form^ to avoid numerical 
integrations that are customary for 
quadrilaterial elements. The approach was also 
used for the closed form evaluation of three 
dimensional hexahedral element integrals^. The 
CPU time used by the closed form evaluation of 
the element integrals has been investigated and 
compared with CPU times required for different 
orders of Gauss numerical integration for 

quadrilateral and hexahedral elements. For the 
quadrilateral elements, the CPU time required to 
evaluate all element integrals is reduced by a 
factor of about 50. For hexahedral elements, 
the time savings are even more significant. 
Although these time savings are beneficial, the 
CPU time required for element integral 
evaluation normally represents less than about 
10% of the total CPU time required for the 
solution. 


Vector Programming Strategies 

The Taylor-Galerkin algorithm was 
implemented with vectorization strategies 
specifically for the Langley VPS 32 (a Cyber 205 
with 32 million words of central memory). This 
computer achieves high computational speed when 
performing operations on long vectors. Vector 
lengths of at least 60 are required to justify 
vectorization efforts with maximum payoff 
achieved for vector lengths of 1000 or more. 
The predominant vector lengths in the 


For an explicit scheme, solutions are 
obtained directly, so that operation (2) 

vectorizes naturally. Operation (3), the 

application of boundary conditions, is an 
intrinsically scalar operation and difficult to 
vectorize. However, use of bit vectors to flag 
boundary nodes and use of VPS 32 FORTRAN 
supplied rountines enables full vectorization of 
this operation. 

A program flow chart for the Taylor- 
Galerkin algorithm is shown in Fig. 2. Note 
that the element Integrals are computed once and 
stored for later use in the transient loop. The 

1 _o 

Lapidus smoothing 1 c represents artificial 
damping used to reduce spurious oscillations. 
The smoothing is typically used only in the 
vicinity of discontinuities, e.g. shock fronts, 
in the solution. 


Application 

Two applications are presented to validate 
and illustrate the basic capabilities of the 
approach. The first application is a one- 
dimensional, transient, linear thermal -stress 
problem for which an exact solution is 

available. The second application is a 
quasistatic, nonlinear, two-dimensional thermal - 
stress problem of a cylinder subject to 
aerodynamic heating. 
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Thermal Stresses in an Elastic Half-Space 

The problem (see Figure 3) consists of an 
elastic half-space x>0 with the plane x=0 free 
of stress, and the medium constrained so that 
the only displacements occur in the x 
direction. The plane x=0 is suddenly at time 
equal zero exposed to convective heating through 
a convection coefficient h. The thermal- 
mechanical coupling term in eq. (lc) is 

neglected to facilitate obtaining an exact 
sol ution. 

The problem was solved using the two-step 
Taylor-Galerkin algorithm with two meshes (Fig. 
3) of quadrilateral elements. Quadri la teral 

elements are not required since the problem is 
one-dimensional, but this application served to 
validate the quadrilateral elements and the 
vector code used for the second application. 
The temperature and stress distributions 
computed for 100 and 500 elements are compared 
with the exact solution in Figs. 4 and 5, 
respectively. 

The computed temperature coincides with the 
exact solution for both meshes. The refined 
mesh was needed to resolve the sharp peak in the 
propagating stress shown in Fig. 5. The coarse 
mesh predicts a smooth stress distribution but 
rounds the sharp peak significantly 
underestimating the peak values. The refined 
mesh solution shows that the solution is 
converging accurately to the exact solution as 
the mesh is refined. 

The problem illustrates the basic 
capability of the algorithm for accurately 
solving transient thermal -stress problems with 
propagating disturbances. Computational times 
for this relatively small problem are modest 
compared to CFD problems. 


Cylinder Subject to Aerodynamics Heating 

The problem (Fig. 6) concerns the thermal- 
structural response of a stainless steel 
cylinder subject to highly-localized aerodynamic 
heating due to supersonic flow. The problem 
simulates the type of aerodynamic heating that 
may occur on leading edges of engine structures 
in hypersonic flight vehicles. The heating is 
time-dependent and results from the interaction 
of two shocks in the flow-field. For 0<t<l s , 
the cylinder is heated symmetrically by a 
supersonic flow that causes a symmetric bow 
shock to form in front of the cylinder leading 
edge. For l<t<2 s , an oblique shock is 

introduced into the flow by rotating the wedge 
shown in Fig. 6. The wedge-induced shock and 
the bow shock then intersect forming a localized 
region of heating at about 10° below the 
cylinder horizontal centerline. Wind tunnel 
experiments suggest that the resulting heating 
can be quite severe. The heating may be 
represented approximately using the convective 
coefficients shown in Fig. 7. The heating 
distribution is modeled using two superimposed 
cosine distributions. The maximum localized 
heating for l<t<2 s is about an order of 

magnitude greater than the heating for 0<t<l s. 


Over the temperature range that occurs in 
the cylinder, the thermal and structural 
properties of the stainless steel cylinder vary 
significantly. The specific heat, thermal 
conductivity and coefficient of thermal 
expansion increase linearly as shown in Fig. 
8. In addition, the stress-strain behavior for 
the material alters significantly with 

temperature as shown in Fig. 9 where it can be 
seen that the elastic moaulus and yield strength 
decrease with increasing temperature. 

The finite element model of the cylinder 
with the thermal and structural boundary 

conditions is shown in Fig. 10. A graded mesh 
of 1500 elements and 1596 nodes is used to model 
one-half of the cylinder. The model represents 
convective heating on the leading edge, and 
assumes negligible heat loss on the other 
cylinder surfaces. The cylinder external 

surfaces are stress-free, but the rigid body 
motion of the cylinder is prohibited by the 
three specified zero displacements shown. Since 
the duration of the heating is much larger than 
the propagation times for thermal -stress waves, 
an excellent approximation is to neglect 

structural inertia effects and perform a 

quasistatic analysis. Thus eqs. (3a) are solved 
by the Taylor-Galerkin algorithm. The 

temperature response is computed in a time- 
accurate manner, but the corresponding stress 
problem is solved independently using the 

fictitious damping constant to march the 

structural response to a steady-state solution 
for a temperature distribution computed at a 
given time. 

The thermal-structural response of the 
cylinder was computed for two cases: (1) 

constant material properties, and (2) 

temperature-dependent properties. Comparative 
temperature distributions for t=ls and t=2s are 
presented in Figs. 11-12, respectively. The 
corresponding stress distributions superimposed 
on greatly exaggerated deformed structures are 
shown in Figs. 13-14. 

The temperatures distribution at 

t=l s (Fig. 11) is changed only slightly by the 
temperature-dependent properties, but the 
temperature distributions at t=2 s (Fig. 12) 
differ considerably. For l<t<2 s, the high 
local heating raises the temperatures near the 
surface significantly. The increasing 

conductivity and specific heat in the 
temperature-dependent properties analysis 
produce smaller gradients and a much lower 
surface temperature. The stress-distribution 
and deformations of Figs. 13 and 14 reflect two 
temperature-dependent property effects. The 
first effect is due to the differences in the 
computed temperatures, and the second is due to 
the difference in the stress-strain behavior of 
the material. Fig. 14, in particular, shows 
these effects dramatical ly. For the constant 
property case. Fig. 14a shows an excessively 
large compressive stress (-325 ksi) well above 
the material's allowable stress. Fig. 14b, 
shows more realistic behavior with smaller 
deformations and stress levels (-150 ksi) still 
high but at acceptable levels. These lower 
levels reflect the reduced temperatures and the 
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inelastic deformations permitted by the 
nonlinear stress-strain curves (Fig. 9). 

This problem reflects some of the 
computational advantages of the explicit Taylor- 
Galerkin algorithm for realistic thermal- 
structural problems. The algorithm is highly 
vectorizable and very large problems can be 
solved on a supercomputer because global 
stiffness matrices are not formed. Non-linear 
property behavior is included conveniently 
through the vectors {E} and {F} independent of 
element stiffness matrices. Element matrices 
are computed from closed-form equations and can 
be computed outside of the transient loop and 
stored for later use. Principal disadvantages 
of the algorithm include its conditional 
stability and tendency to produce oscillatory 
solutions in the presence of sharp solution 
discontinuities. These disadvantages have been 
encountered in inviscid flow computations 1 * 2 and 
nonlinear thermal -structural applications with 
propagating shock waves may show similar 
results. Further computations are needed to 
investigate these possibilities. Another area 
worthy of further study is methods for 
accelerating the convergence of the time- 
marching solution to the steady-state structural 
equilibrium problem associated with quasistatic 
thermal-structural analysis. 


Concluding Remarks 

A Taylor-Galerkin finite element solution 
algorithm for transient nonlinear thermal- 
structural analysis of large, complex structural 
problems Is described. The two-step Taylor- 
Galerkin algorithm is an application of an 
algorithm recently developed for problems in 
compressible fluid dynamics. Two thermal- 

structural formulations are described. The 
first formulation is for thermal-structural 
problems where a mechanical coupling term is 
present in the heat transfer energy equation and 
the Inertial terms are retained in the 
structural equations of motion. The second 
formulation is for a quasistatic problem where 
thermal-mechanical coupling and structural 
inertia terms are neglected. Solutions to 

transient and equilibrium problems are obtained 
by time-marching. For a lumped mass approach, 
the algorithm is conditionally stable. The 
algorithm has been implemented on the VPS-32 
vector computer with special programming 
strategies to yield very high computational 
speeds. 

Two applications of the algorithm to 
thermal -structural problems are presented. The 
first application is a one-dimensional, 
transient linear thermal -stress problem for an 
elastic half-space. Comparisons of finite 
element predictions with an exact solution 
validate the approach and illustrate the 
capability to capture propagating stress waves 
accurately. The second application is a 
quasistatic, nonlinear two-dimensional thermal- 
stress problem of a cylinder subject to high 
localized aerodynamic heating. Comparisons of 
constant material property and variable property 


solutions illustrated the importance of 
nonlinear effects in realistic problems and the 
capability of the algorithm to incorporate the 
nonlinearities effectively. 

The applications have validated the 
fundamental capabilities of the algorithm for 
two basic thermal-structural formulations. 
Additional study is needed for more demanding 
transient, nonlinear thermal-stress problems. 
Further experience is needed also for 
quasistatic problems particularly with methods 
to accelerate the convergence of the steady- 
state structural problem. Since the algorithm 
can be applied to fluid, thermal and structural 
problems there is potential for developing, with 
one methodology, the capability fpr integrating 
these analyses into a single program. Such an 
integrated methodology will permit the first 
computational solution to problems with strong 
fluid-thermal-structural interactions. The 
development of an integrated fluid- thermal- 
structural analysis capability is a current 
research goal. 
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STRUCTURAL MODEL 


Fig. 1 Two dimensional finite element thermal- 
structural model. 


BOUNDARY CONDITIONS 

• Convective heating 

k 10. tl = hl TIO.II - T^l 

• Zero surface traction 
0 , 10.11 =0 



INITIAL CONDITIONS 

Tlx. 01 = Tg I T.p 

u lx. 01 = 0 T 

u (x 0] =0 ' 1 


Fig. 3 Transient thermal stress in an elastic 
half space. 



X 


Fig. 4 Comparative temperature distribution at 
time t * 5. 



Fig. 2 Program flowchart for Taylor-Galerkin 
algorithm. 



X 


Fig. 5 Comparative thermal stress distributions 
at time t = 5. 







ORIGINAL PAGE IS 
OF POOR QUALITY 



Fig. 6 Thermal-structural analysis of a Fig. 9 Temperature dependent stress-strain 

cylinder under shock wave heating. curves of cylinder. 



Fig. 7 Convective coefficient distributions on 
cylindrical surface. 


Fig. 10 Finite element thermal-structural model 
for cylinder. 



T. R 


Fig. 11 Comparative temperature distributions 
(in R) at 1 second. 


Fig. 8 Temperature dependent material 
properties of cylinder. 
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(a) Constant material properties. (b) Temperature dependent 

material properties. 

Fig. 12 Comparative temperature distributions 
(in R) at 2 seconds. 




(a) Constant material properties, (b) Temperature dependent material 
Elastic analysis. properties. Inelastic analysis. 


Fig. 13 Comparative c 
distributions 



(a) Constant material properties. 
Elastic analysis. 


rcumferential stress 
(in ksi) at 1 second. 



(b) Temperature dependent material 
properties. Inelastic analysis. 


Fig. 14 Comparative circumferential stress 
distributions (in ksi) at 2 seconds. 
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